{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 在量子化学计算中应用量子变分求解器\n",
    "\n",
    "[![下载Notebook](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/master/resource/_static/logo_notebook.png)](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/notebook/master/mindquantum/zh_cn/mindspore_vqe_for_quantum_chemistry.ipynb)&emsp;[![下载样例代码](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/master/resource/_static/logo_download_code.png)](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/notebook/master/mindquantum/zh_cn/mindspore_vqe_for_quantum_chemistry.py)&emsp;[![查看源文件](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/master/resource/_static/logo_source.png)](https://gitee.com/mindspore/docs/blob/master/docs/mindquantum/docs/source_zh_cn/vqe_for_quantum_chemistry.ipynb)\n",
    "\n",
    "## 概述\n",
    "\n",
    "量子化学，指的是运用量子力学的基本理论及方法，求解含时或定态薛定谔方程的数值解。在高性能计算机上进行量子化学模拟已成为研究材料的物理、化学性质的重要手段。然而，精确求解薛定谔方程具有指数级的复杂度，可模拟的化学体系规模严重受制于此。近年量子计算的发展为解决这个问题提供了一条可行的路，有望在量子计算机上实现多项式复杂度下对薛定谔方程的高精度求解。\n",
    "\n",
    "[Peruzzo等人](https://doi.org/10.1038/ncomms5213)在2014年首次将量子变分求解器(Variational quantum eigensolver, VQE)结合[幺正耦合簇理论](https://linkinghub.elsevier.com/retrieve/pii/S0009261489873725)用于量子化学的模拟中，实现了He-H<sup>+</sup>基态能量的求解。量子变分求解器是一个量子--经典混合算法，在基于量子算法的化学模拟中应用广泛，本教程将介绍使用量子变分求解器求解分子体系基态能量的方法。\n",
    "\n",
    "本教程的主要内容包括如下几个部分：\n",
    "\n",
    "1. 量子化学原理简介。\n",
    "2. 量子变分求解器的应用。\n",
    "3. 使用MindQuantum实现高效自动求导的VQE模拟。\n",
    "\n",
    "> 本文档适用于CPU环境。\n",
    ">\n",
    "> 你可以在这里找到完整的可运行的样例代码：<https://gitee.com/mindspore/mindquantum/blob/master/tutorials/source/7.vqe_for_quantum_chemistry.py>。\n",
    "\n",
    "## 环境准备\n",
    "\n",
    "本教程需要安装以下环境：\n",
    "\n",
    "- NumPy\n",
    "- SciPy\n",
    "- [mindquantum](https://gitee.com/mindspore/mindquantum)\n",
    "- [mindspore](https://gitee.com/mindspore/mindspore)\n",
    "- PySCF\n",
    "- openfermion\n",
    "- openfermionpyscf\n",
    "\n",
    "> 以上依赖都可通过`pip`命令来安装。\n",
    "\n",
    "## 导入依赖\n",
    "\n",
    "导入本教程所依赖模块"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from openfermion.chem import MolecularData\n",
    "from openfermionpyscf import run_pyscf\n",
    "from mindquantum import Circuit, X, Hamiltonian, Simulator\n",
    "from mindquantum.algorithm import generate_uccsd\n",
    "import mindspore as ms\n",
    "import mindspore.context as context\n",
    "from mindspore import Parameter\n",
    "\n",
    "context.set_context(mode=context.PYNATIVE_MODE, device_target=\"CPU\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 量子化学计算方法\n",
    "\n",
    "量子化学的核心问题在于求解薛定谔方程（Schrödinger Equation）。一般来说，求解含时薛定谔方程（Time-dependent Schrödinger Equation)较为复杂，故引入玻恩-奥本海默近似（Born-Oppenheimer approximation, BO approximation）。BO近似认为，原子核质量远大于电子、运动速度远低于电子，故可以将两者进行分离变量，单独讨论原子核或电子的运动，于是可得到如下不含时的电子运动方程，也称为定态薛定谔方程：\n",
    "\n",
    "$$\n",
    "\\hat{H} |\\Psi\\rangle = E |\\Psi\\rangle\n",
    "$$\n",
    "\n",
    "其中$\\hat{H}$包含以下三项：\n",
    "\n",
    "$$\n",
    "\\hat{H} = \\hat{K} _{e} + \\hat{V} _{ee} + \\hat{V} _{Ne}\n",
    "$$\n",
    "\n",
    "分别为电子动能、电子-电子势能和电子-核势能。\n",
    "\n",
    "有多种数值算法可以求解定态薛定谔方程。本教程将介绍其中的一类：波函数方法。波函数方法直接求解给定分子哈密顿量的本征波函数和本征能量，目前有大量的开源软件包可实现，如[PySCF](http://pyscf.org/)等。此处从一个简单的例子：氢化锂分子开始，使用openfermion结合openfermionpyscf插件进行。首先定义分子结构："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geometry: \n",
      " [['Li', [0.0, 0.0, 0.0]], ['H', [0.0, 0.0, 1.5]]]\n"
     ]
    }
   ],
   "source": [
    "dist = 1.5\n",
    "geometry = [\n",
    "    [\"Li\", [0.0, 0.0, 0.0 * dist]],\n",
    "    [\"H\", [0.0, 0.0, 1.0 * dist]],\n",
    "]\n",
    "basis = \"sto3g\"\n",
    "spin = 0\n",
    "print(\"Geometry: \\n\", geometry)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上面的代码定义了一个Li-H键长为1.5Å分子。使用STO-3G基组进行计算。接下来使用openfermionpyscf，调用PySCF进行HF、CCSD和FCI计算。这三种方法属于波函数方法，开始计算之前，先对这些方法作一个简单的介绍。\n",
    "\n",
    "### 波函数方法\n",
    "\n",
    "求解定态薛定谔方程的方法之一是[Hartree-Fock（HF）](https://doi.org/10.1098/rspa.1935.0085)方法，该方法在二十世纪三十年代左右由Hartree等人提出，是量子化学计算中的基本方法。HF方法引入了单行列式近似，即$N$-电子体系的波函数由一个行列式形式的波函数表示：\n",
    "\n",
    "$$\n",
    "| \\Psi \\rangle = | \\psi_{1} \\psi_{2} \\psi_{3} \\dots \\psi_{N} \\rangle\n",
    "$$\n",
    "\n",
    "其中$| \\psi_{1} \\psi_{2} \\psi_{3} \\dots \\rangle$代表由一组自旋轨道波函数$\\{ \\pi_{i} \\}$构成的N阶行列式。\n",
    "自旋轨道波函数$\\psi_{i}$可进一步用一组形式已知的基函数展开：\n",
    "\n",
    "$$\\psi_{i} = \\phi_{i} \\eta_{i}$$\n",
    "\n",
    "$$\\phi_{i} = \\sum_{\\mu}{C_{\\mu i} \\chi_{\\mu}}$$\n",
    "\n",
    "其中$\\{\\chi_{\\mu}\\}$被称为基函数，可以是高斯函数等。\n",
    "该近似考虑了电子间的交换作用，但是忽略了电子间的关联作用，故无法正确计算如解离能等性质。\n",
    "\n",
    "HF方法的改进可以从波函数展开定理出发。波函数展开定理可以表述为，若$\\{ \\psi_{i} \\}$是一组完备的自旋轨道波函数，则$N$-电子体系波函数可以由$\\{ \\psi_{i} \\}$构成的行列式波函数精确展开：\n",
    "\n",
    "$$\n",
    "| \\Psi \\rangle = \\sum^{\\infty} _ {i_{1} < i_{2} < \\dots < i_{N}} {C_{i_{1} i_{2} \\dots i_{N}} | \\psi_{i_{1}} \\psi_{i_{2}} \\dots \\psi_{i_{N}} \\rangle}\n",
    "$$\n",
    "\n",
    "由此可得到Configuration Interaction（CI）方法：\n",
    "\n",
    "$$\n",
    "| \\Psi_{CI} \\rangle = C_{0} | \\Psi_{HF} \\rangle + \\sum^{a\\rightarrow\\infty} _{i\\in occ\\\\\\\\a\\not\\in occ}{C^{a} _{i} | \\Psi^{a} _{i} \\rangle } + \\sum^{ab\\rightarrow\\infty} _{ij\\in occ\\\\\\\\ab\\not\\in occ}{C^{ab} _{ij} | \\Psi^{ab} _{ij} \\rangle }\n",
    "$$\n",
    "\n",
    "上式中的$| \\Psi^{a}_{i} \\rangle + \\dots$代表电子由轨道$i$激发到轨道$a$的单激发波函数，以此类推。只考虑单激发和双激发的CI被称为CISD，即Configuration Interaction with singles and doubles。将基态HF波函数一直到N激发波函数全部考虑在内的Configuration Interaction被称为Full Configuration Interaction（FCI），FCI波函数是定态薛定谔方程在给定基函数下的精确解。\n",
    "\n",
    "### 二次量子化\n",
    "\n",
    "在二次量子化表述下，体系的哈密顿量具有如下形式：\n",
    "\n",
    "$$\n",
    "\\hat{H} = \\sum_{p, q}{h^{p} _ {q} E^{p} _ {q}} + \\sum_{p, q, r, s}{\\frac{1}{2} g^{pq} _ {rs} E^{pq} _ {rs} }\n",
    "$$\n",
    "\n",
    "其中$E^{p}_{q}$和$E^{pq}_{rs}$分别为：\n",
    "\n",
    "$$\n",
    "E^{p} _ {q} = a^{\\dagger} _ {p} a_{q}\n",
    "\n",
    "$$\n",
    "\n",
    "$$\n",
    "E^{pq} _ {rs} = a^{\\dagger} _ {p} a^{\\dagger} _ {q} a_{r} a_{s}\n",
    "$$\n",
    "\n",
    "$a^{\\dagger}_{p}$和$a_{q}$分别为产生算符（Creation Operator）和湮灭算符（Annihilation Operator）。\n",
    "\n",
    "使用二次量子化的表述方法，可以非常方便地表示激发态波函数：\n",
    "\n",
    "$$\n",
    "| \\Psi^{abc\\dots} _ {ijk\\dots} \\rangle = a^{\\dagger} _ {a} a^{\\dagger} _ {b} a^{\\dagger} _ {c} \\dots a_{i} a_{j} a_{k} \\dots | \\Psi \\rangle\n",
    "$$\n",
    "\n",
    "CI方法的一个改进是耦合簇理论（Coupled-Cluster theory, CC）。CC引入指数化算符：\n",
    "\n",
    "$$\n",
    "| \\Psi_{CC} \\rangle = \\exp{(\\hat{T})} | \\Psi_{HF} \\rangle\n",
    "$$\n",
    "\n",
    "其中耦合簇算符$\\hat{T}$为对激发算符的求和：\n",
    "\n",
    "$$\n",
    "\\hat{T} = \\sum_{p\\not\\in occ\\\\\\\\q\\in occ}{\\theta^{p} _ {q} E^{p} _ {q}} + \\sum_{pq\\not\\in occ\\\\\\\\rs\\in occ}{\\theta^{pq} _ {rs} E^{pq} _ {rs}} + \\dots\n",
    "$$\n",
    "\n",
    "其中$\\theta$和CI方法中的$C$类似，是待求解的参数。由指数的泰勒展开易知，即使耦合簇算符$\\hat{T}$中只包含低阶激发项，$\\exp{(\\hat{T})}$也可以隐含部分高阶激发，这也使得CC方法向FCI波函数收敛的速度要远快于CI，同样截断到K激发，如K=2，CCSD的精度会超过CISD。\n",
    "\n",
    "<!--\n",
    "一般而言，若一个方法可以达到化学精度，即由此方法计算的能量和FCI能量之间的差值小于1 kcal/mol，则认为这个方法具有良好的精度，截断到三激发的CCSD(T)在大部分情况下都能符合这个标准\n",
    "-->\n",
    "\n",
    "电子关联作用的效果是使得总能量降低，故HF得到的基态能量会略高于CCSD和FCI。另外，从上述理论不难发现，FCI的计算量远大于CCSD和HF。我们使用openfermion封装的`MolecularData`和openfermionpyscf封装的`run_pyscf`函数来进行演示："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hartree-Fock energy:  -7.8633576215351164 Ha\n",
      "CCSD energy:  -7.8823529091527007 Ha\n",
      "FCI energy:  -7.8823622867987213 Ha\n"
     ]
    }
   ],
   "source": [
    "molecule_of = MolecularData(\n",
    "    geometry,\n",
    "    basis,\n",
    "    multiplicity=2 * spin + 1\n",
    ")\n",
    "molecule_of = run_pyscf(\n",
    "    molecule_of,\n",
    "    run_scf=1,\n",
    "    run_ccsd=1,\n",
    "    run_fci=1\n",
    ")\n",
    "\n",
    "print(\"Hartree-Fock energy: %20.16f Ha\" % (molecule_of.hf_energy))\n",
    "print(\"CCSD energy: %20.16f Ha\" % (molecule_of.ccsd_energy))\n",
    "print(\"FCI energy: %20.16f Ha\" % (molecule_of.fci_energy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在上面的例子中，我们运行了Hartree-Fock（HF）、CCSD、FCI进行总能量的计算。若对运行时间进行统计，会发现$T_{HF}<T_{CCSD}\\ll T_{FCI}$，换成计算量更大的体系如乙烯分子等会更明显一些。此外，对于计算得到的总能量，有$E_{HF}>E_{CCSD}>E_{FCI}$。计算完成后，我们将结果保存到`molecule_file`文件（即`molecule_of.filename`）中："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/xuxs/anaconda3/lib/python3.8/site-packages/openfermion/testing/data/H1-Li1_sto3g_singlet\n"
     ]
    }
   ],
   "source": [
    "molecule_of.save()\n",
    "molecule_file = molecule_of.filename\n",
    "print(molecule_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "量子化学计算的一大阻碍是计算量。随着体系大小（电子数、原子数）的增加，求解FCI波函数和基态能量的时间消耗大约以$2^{N}$增长，即使是较小的分子如乙烯分子等，进行FCI计算也并不容易。量子计算机的出现为此提供了一条可能的解决途径，已有的研究表明，量子计算机可以多项式的时间复杂度模拟哈密顿量的含时演化，在量子处理器上进行化学模拟相较于经典计算机有指数级的加速。本教程将介绍其中一类量子算法：量子变分求解器。\n",
    "\n",
    "## 量子变分求解器\n",
    "\n",
    "量子变分求解器（Variational Quantum Eigensolver, VQE）是一类量子-经典混合（Hybrid quantum-classical）算法，应用变分原理实现对基态波函数的求解。其中，变分参数的优化步在经典计算机上进行。\n",
    "\n",
    "### 变分原理\n",
    "\n",
    "变分原理可使用如下形式表述：\n",
    "\n",
    "$$\n",
    "E_{0} \\le \\frac{\\langle \\Psi_{t} | \\hat{H} | \\Psi_{t} \\rangle}{\\langle \\Psi_{t} | \\Psi_{t} \\rangle}\n",
    "$$\n",
    "\n",
    "上式中的$| \\Psi_{t} \\rangle$代表试探波函数。变分原理表明，在满足一定的条件下，任意试探波函数得到的基态能量总是大于等于真实的基态能量。变分原理为求解分子基态薛定谔方程提供了一种方法：使用一个参数化的函数$f(\\theta)$作为精确基态波函数的近似，通过优化参数$\\theta$来逼近精确的基态能量。\n",
    "\n",
    "### 初态制备\n",
    "\n",
    "在二次量子化表述下，$N$-电子HF波函数也具有非常简洁的形式：\n",
    "\n",
    "$$\n",
    "| \\Psi_{HF} \\rangle = \\prod^{i=0} _{N-1}{a^{\\dagger} _{i}| 0 \\rangle}\n",
    "$$\n",
    "\n",
    "上式搭建了一个由量子化学波函数到量子计算的桥梁：用$|0\\rangle$代表非占据轨道，用$|1\\rangle$代表电子占据的轨道，由此可以将$N$-电子HF波函数映射为由一串$M+N$个量子比特$| 00\\dots 11\\dots \\rangle$，$M$代表非占据轨道的数量。\n",
    "\n",
    "以下代码构造了对应于LiH分子的HF初态波函数。在Jordan-Wigner变换下，相当于将$N$个$\\text{X}$门作用于$|000\\dots\\rangle$上。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "q0: ──X──\n",
      "         \n",
      "q1: ──X──\n",
      "         \n",
      "q2: ──X──\n",
      "         \n",
      "q3: ──X──\n"
     ]
    }
   ],
   "source": [
    "hartreefock_wfn_circuit = Circuit([X.on(i) for i in range(molecule_of.n_electrons)])\n",
    "print(hartreefock_wfn_circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "基于此，我们可以构造如下形式的试探波函数：\n",
    "\n",
    "$$\n",
    "| \\Psi_{t} \\rangle = U(\\theta) | \\Psi_{HF} \\rangle\n",
    "$$\n",
    "\n",
    "其中$U(\\theta)$代表一个可通过量子线路模拟的幺正变换，$| \\Psi_{HF} \\rangle$作为初态，可通过多个单比特$\\text{X}$门来方便地制备。$U(\\theta) | \\Psi_{HF} \\rangle$的具体形式也被称为波函数拟设。\n",
    "\n",
    "### 波函数拟设\n",
    "\n",
    "前文提到的耦合簇理论是一个非常高效的波函数拟设。在量子计算机上使用，需要作一些修改：\n",
    "\n",
    "$$\n",
    "| \\Psi_{UCC} \\rangle = \\exp{(\\hat{T} - \\hat{T}^{\\dagger})} | \\Psi_{HF} \\rangle\n",
    "$$\n",
    "\n",
    "UCC即幺正耦合簇(Unitary Coupled-Cluster theory)，$\\hat{T}^{\\dagger}$代表$\\hat{T}$的厄米共轭。如此，$\\exp{(\\hat{T} - \\hat{T}^{\\dagger})}$即为幺正算符。[Peruzzo等人](https://doi.org/10.1038/ncomms5213)在2014年首次使用VQE结合UCCSD(Unitary coupled-cluster with singles and doubles)拟设进行了量子计算机上的化学模拟实验。值得注意的是幺正耦合簇默认了耦合簇算符中的参数$\\{\\theta\\}$是实数。在分子体系中该假设不会有问题；在周期性体系中，[刘杰等人](https://doi.org/10.1021/acs.jctc.0c00881)的研究表明幺正耦合簇会因为忽略复数部分而造成误差。本教程暂时不讨论幺正耦合簇在周期性体系中的应用。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用mindquantum的circuit模块中的`generate_uccsd`函数可读取先前保存在`molecule_file`的计算结果，“一键”构造UCCSD波函数拟设，以及其对应的量子线路："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ccsd:-7.882352909152701.\n",
      "fci:-7.882362286798721.\n"
     ]
    }
   ],
   "source": [
    "ansatz_circuit, \\\n",
    "init_amplitudes, \\\n",
    "ansatz_parameter_names, \\\n",
    "hamiltonian_QubitOp, \\\n",
    "n_qubits, n_electrons = generate_uccsd(molecule_file, th=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`generate_uccsd`将幺正耦合簇相关的函数打包了起来，包括导出分子哈密度量、构造幺正耦合簇拟设算符、提取CCSD计算的耦合簇系数等多个步骤。该函数通过输入分子的文件路径来读取该分子，参数`th`是表示量子线路中哪些参数需要更新梯度的阈值。在[分步构造幺正耦合簇拟设](#step-by-step)章节，我们会演示如何使用mindquantum的相关接口分步完成其中包含的步骤。完整的量子线路包含HF初态+UCCSD拟设，如下代码所示："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "============================Circuit Summary============================\n",
      "|Total number of gates  : 15172.                                      |\n",
      "|Parameter gates        : 640.                                        |\n",
      "|with 44 parameters are : p0, p8, p1, p9, p2, p10, p3, p11, p4, p12...|\n",
      "|Number qubit of circuit: 12                                          |\n",
      "=======================================================================\n",
      "Number of parameters: 44\n"
     ]
    }
   ],
   "source": [
    "total_circuit = hartreefock_wfn_circuit + ansatz_circuit\n",
    "total_circuit.summary()\n",
    "print(\"Number of parameters: %d\" % (len(ansatz_parameter_names)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于LiH分子而言，其UCCSD波函数拟设中包含44个变分参数。该线路总共的量子比特门数量为12612，总共需要12个量子比特进行模拟。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### VQE的一般流程\n",
    "\n",
    "使用VQE进行分子基态求解的一般流程如下：\n",
    "\n",
    "1. 制备HF初态：$| 00\\dots11\\dots \\rangle$；\n",
    "2. 定义波函数拟设，如UCCSD等；\n",
    "3. 将波函数拟设转化为参数化的量子线路；\n",
    "4. 初始化变分参数，如全设为0等；\n",
    "5. 在量子计算机上多次测量得到分子哈密顿量在该套变分参数下的能量$E(\\theta)$以及能量关于参数的导数$\\{ {\\partial E} / {\\partial \\theta_{i}} \\}$\n",
    "6. 在经典计算机上使用优化算法，如梯度下降、BFGS等更新变分参数；\n",
    "7. 将新的变分参数传入量子线路中进行更新；\n",
    "8. 重复步骤(5)到(7)，直到满足收敛标准；\n",
    "9. 结束\n",
    "\n",
    "在第5步中，求取能量关于参数的导数$\\{ {\\partial E} / {\\partial \\theta_{i}} \\}$在量子计算机上可通过parameter-shift rule来进行，在模拟器中也可通过模拟parameter-shift rule或者有限差分法来计算，是个较为耗时的过程。mindquantum基于mindspore框架，提供了类似于机器学习的自动求导功能，可以在模拟中可以高效计算变分量子线路的导数。以下使用mindquantum构造带自动求导功能的参数化UCCSD量子线路："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = Simulator('projectq', total_circuit.n_qubits)\n",
    "molecule_pqc = sim.get_expectation_with_grad(Hamiltonian(hamiltonian_QubitOp), total_circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通过将参数的具体数值传入`molecule_pqc`，即可得到对应于此变分参数的能量$E(\\theta)=\\langle \\Psi_{UCC}(\\theta) | \\hat{H} | \\Psi_{UCC}(\\theta) \\rangle$以及关于每个变分参数的导数。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来需要进行VQE优化的(5)~(7)步，即对变分量子线路进行优化。我们可以借助MindQuantum框架，使用变分量子线路算子`molecule_pqc`构造一个神经网络模型，然后通过类似于训练神经网络的方法来优化变分参数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mindquantum.framework import MQAnsatzOnlyLayer\n",
    "\n",
    "molecule_pqcnet = MQAnsatzOnlyLayer(molecule_pqc, 'Zeros')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "此处我们构造了一个基本的`MQAnsatzOnlyLayer`作为模型示例，该模型可以和常规的机器学习模型类似使用，比如优化权重、计算导数等"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "构造的`MQAnsatzOnlyLayer`使用`\"Zeros\"`关键字，将所有的变分参数初始化为0。使用CCSD（耦合簇理论）或者MP2（二阶多体微扰论）的计算结果也可以作为幺正耦合簇变分参数的初始值。此时有$E(\\vec{0})=\\langle \\Psi_{UCC}(\\vec{0}) | \\hat{H} | \\Psi_{UCC}(\\vec{0}) \\rangle = E_{HF}$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial energy:  -7.8633575439453125\n"
     ]
    }
   ],
   "source": [
    "initial_energy = molecule_pqcnet()\n",
    "print(\"Initial energy: %20.16f\" % (initial_energy.asnumpy()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后使用mindspore的Adam优化器进行优化，学习率设置为$1\\times 10^{-2}$，优化终止标准设置为$\\left.|\\epsilon|\\right. = \\left.|E^{k+1} - E^{k}|\\right. \\le 1\\times 10^{-8}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Step   0 energy  -7.8633575439453125\n",
      "Step   5 energy  -7.8726239204406738\n",
      "Step  10 energy  -7.8821778297424316\n",
      "Step  15 energy  -7.8822836875915527\n",
      "Step  20 energy  -7.8823199272155762\n",
      "Step  25 energy  -7.8823370933532715\n",
      "Step  30 energy  -7.8815641403198242\n",
      "Step  35 energy  -7.8786268234252930\n",
      "Step  40 energy  -7.8778734207153320\n",
      "Step  45 energy  -7.8808088302612305\n",
      "Step  50 energy  -7.8819031715393066\n",
      "Step  55 energy  -7.8821558952331543\n",
      "Step  60 energy  -7.8823504447937012\n",
      "Optimization completed at step  63\n",
      "Optimized energy:  -7.8823523521423340\n",
      "Optimized amplitudes: \n",
      " [ 2.40339446e-04  1.89154677e-03  3.49554531e-02  1.59917790e-02\n",
      "  2.33248898e-07  9.09393420e-04 -1.79268172e-05  1.41595434e-02\n",
      "  6.28582342e-08  9.08669957e-04 -1.49387897e-05  1.41652254e-02\n",
      " -5.46666037e-04  4.26779327e-04  2.86067789e-03  5.38198128e-02\n",
      "  2.32545775e-04 -2.78862785e-07 -7.10907813e-08 -7.98562283e-08\n",
      "  7.00364581e-07 -9.21200325e-08 -6.73263187e-08  1.26236855e-05\n",
      " -1.04519488e-04  7.97090179e-04 -4.01437364e-06 -3.34858555e-06\n",
      " -5.49289174e-02  3.09006264e-03  7.01365061e-05 -1.36400865e-06\n",
      " -1.35536197e-06  4.63907739e-08  5.32547162e-08 -2.34681625e-08\n",
      "  3.92657455e-07  5.11744884e-06  3.09006032e-03 -2.05122589e-07\n",
      "  5.91138871e-08 -2.44064164e-08  4.26194856e-06  3.72134935e-04]\n"
     ]
    }
   ],
   "source": [
    "optimizer = ms.nn.Adagrad(molecule_pqcnet.trainable_params(), learning_rate=4e-2)\n",
    "train_pqcnet = ms.nn.TrainOneStepCell(molecule_pqcnet, optimizer)\n",
    "\n",
    "eps = 1.e-8\n",
    "energy_diff = eps * 1000\n",
    "energy_last = initial_energy.asnumpy() + energy_diff\n",
    "iter_idx = 0\n",
    "while abs(energy_diff) > eps:\n",
    "    energy_i = train_pqcnet().asnumpy()\n",
    "    if iter_idx % 5 == 0:\n",
    "        print(\"Step %3d energy %20.16f\" % (iter_idx, float(energy_i)))\n",
    "    energy_diff = energy_last - energy_i\n",
    "    energy_last = energy_i\n",
    "    iter_idx += 1\n",
    "\n",
    "print(\"Optimization completed at step %3d\" % (iter_idx - 1))\n",
    "print(\"Optimized energy: %20.16f\" % (energy_i))\n",
    "print(\"Optimized amplitudes: \\n\", molecule_pqcnet.weight.asnumpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以看到，幺正耦合簇给出的计算结果和FCI非常接近，具有良好的精度。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 分步构造幺正耦合簇拟设\n",
    "\n",
    "<a id=\"step-by-step\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在上文中，我们使用了`generate_uccsd`一步构造出了幺正耦合簇拟设所需要的所有内容，此处我们将步骤拆分，分别得到我们需要的耦合簇算符、对应的量子线路以及取自于经典CCSD计算结果的变分参数初猜值。\n",
    "首先，导入部分额外依赖，主要包含mindquantum中hiqfermion模块的相关函数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mindquantum.algorithm.nisq.chem import Transform\n",
    "from mindquantum.algorithm.nisq.chem import get_qubit_hamiltonian\n",
    "from mindquantum.algorithm.nisq.chem import uccsd_singlet_generator, uccsd_singlet_get_packed_amplitudes\n",
    "from mindquantum.core.operators import TimeEvolution\n",
    "from mindquantum.framework import MQAnsatzOnlyLayer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "分子哈密顿量使用`get_qubit_hamiltonian`，读取之前的计算结果得到："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "hamiltonian_QubitOp = get_qubit_hamiltonian(molecule_of)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于幺正耦合簇算符$ \\hat{T} - \\hat{T}^{\\dagger} $，可以使用`uccsd_singlet_generator`进行构造。提供总量子比特数（总自旋轨道数）和总电子数，并设置参数`anti_hermitian=True`："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "ucc_fermion_ops = uccsd_singlet_generator(\n",
    "    molecule_of.n_qubits, molecule_of.n_electrons, anti_hermitian=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上一步构造的`ucc_fermion_ops`是参数化的。使用Jordan-Wigner变换将费米子激发算符映射为Pauli算符："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "ucc_qubit_ops = Transform(ucc_fermion_ops).jordan_wigner()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来，我们需要得到幺正算符 $ \\exp{(\\hat{T} - \\hat{T}^{\\dagger})} $ 所对应的量子线路。`TimeEvolution`可生成$ \\exp{(-i\\hat{H}t)} $所对应的线路，其中$ \\hat{H} $是一个厄米算符，$t$是实数。需要注意的是，使用`TimeEvolution`时，`ucc_qubit_ops`中已经包含了复数因子$i$，所以我们需要将`ucc_qubit_ops`除以$i$，或者提取其虚部："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "ansatz_circuit = TimeEvolution(ucc_qubit_ops.imag, 1.0).circuit\n",
    "ansatz_parameter_names = ansatz_circuit.params_name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们使用`ansatz_parameter_names`记录该线路中的参数名。到目前为止，我们已经得到了VQE量子线路所需要内容，包括哈密顿量`hamiltonian_QubitOp`、参数化的波函数拟设线路`ansatz_circuit`，故可仿照前文，得到完整的态制备线路。其中Hartree-Fock参考态复用之前的`hartreefock_wfn_circuit`："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==================================Circuit Summary==================================\n",
      "|Total number of gates  : 15172.                                                  |\n",
      "|Parameter gates        : 640.                                                    |\n",
      "|with 44 parameters are : s_0, d1_0, s_1, d1_1, s_2, d1_2, s_3, d1_3, s_4, d1_4...|\n",
      "|Number qubit of circuit: 12                                                      |\n",
      "===================================================================================\n"
     ]
    }
   ],
   "source": [
    "total_circuit = hartreefock_wfn_circuit + ansatz_circuit\n",
    "total_circuit.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下一步，需要为变分参数提供一个合理的初始值。前文构造的`PQCNet`默认使用0作为初猜，在大多数情况下是可行的。不过，使用CCSD的计算数据作为UCC的出发点，可能会有更好的结果。使用`uccsd_singlet_get_packed_amplitudes`函数从`molecule_of`提取CCSD的参数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_amplitudes_ccsd = uccsd_singlet_get_packed_amplitudes(\n",
    "    molecule_of.ccsd_single_amps, molecule_of.ccsd_double_amps, molecule_of.n_qubits, molecule_of.n_electrons)\n",
    "init_amplitudes_ccsd = [init_amplitudes_ccsd[param_i] for param_i in ansatz_parameter_names]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用`MQAnsatzOnlyLayer`可以方便地由参数、量子线路获得以变分量子线路为基础的机器学习模型："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "grad_ops = Simulator('projectq', total_circuit.n_qubits).get_expectation_with_grad(\n",
    "    Hamiltonian(hamiltonian_QubitOp.real),\n",
    "    total_circuit)\n",
    "\n",
    "molecule_pqcnet = MQAnsatzOnlyLayer(grad_ops)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用`init_amplitudes_ccsd`（即CCSD计算的耦合簇系数）作为初始变分参数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial energy:  -7.8173098564147949\n"
     ]
    }
   ],
   "source": [
    "molecule_pqcnet.weight = Parameter(ms.Tensor(init_amplitudes_ccsd, molecule_pqcnet.weight.dtype))\n",
    "initial_energy = molecule_pqcnet()\n",
    "print(\"Initial energy: %20.16f\" % (initial_energy.asnumpy()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这个例子中，CCSD初猜并没有提供一个更好的起点。读者可以对更多的分子、更多种类的初始值（如随机数初猜）等进行测试和探究。最后进行VQE的优化步骤，优化器依然使用Adam，收敛标准不变。优化所用的代码与前文基本一致，注意更新相应的变量即可："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eps:  1e-08\n",
      "Step   0 energy  -7.8173098564147949\n",
      "Step   5 energy  -7.8740758895874023\n",
      "Step  10 energy  -7.8818783760070801\n",
      "Step  15 energy  -7.8821649551391602\n",
      "Step  20 energy  -7.8822622299194336\n",
      "Step  25 energy  -7.8823080062866211\n",
      "Step  30 energy  -7.8822288513183594\n",
      "Step  35 energy  -7.8758554458618164\n",
      "Step  40 energy  -7.8761253356933594\n",
      "Step  45 energy  -7.8807921409606934\n",
      "Step  50 energy  -7.8818383216857910\n",
      "Step  55 energy  -7.8821811676025391\n",
      "Step  60 energy  -7.8823504447937012\n",
      "Optimization completed at step  63\n",
      "Optimized energy:  -7.8823523521423340\n",
      "Optimized amplitudes: \n",
      " [-2.42472161e-04  1.89258391e-03 -3.46013680e-02  1.59353409e-02\n",
      " -8.40432079e-08  9.09362687e-04  2.01798011e-05  1.41534032e-02\n",
      " -2.81526667e-07  9.08639806e-04  2.00776722e-05  1.41590768e-02\n",
      "  5.45396935e-04  4.26715094e-04 -2.84755090e-03  5.38354851e-02\n",
      "  2.29954778e-04  9.55212727e-07 -1.24844689e-07 -9.20767249e-08\n",
      " -4.53033465e-07 -7.44455733e-08 -8.83169875e-08  1.17984437e-05\n",
      " -1.04996754e-04  7.94677646e-04 -4.50417019e-06 -4.49753043e-06\n",
      " -5.48430867e-02  3.08870710e-03 -6.50319926e-05  1.26427835e-06\n",
      "  1.25660222e-06 -2.82077963e-07  7.96948143e-08 -3.28978906e-08\n",
      " -3.63568660e-07  5.76087541e-06  3.08870478e-03  8.82309266e-08\n",
      "  5.73797401e-08 -2.53652850e-08  5.72846511e-06  3.71275470e-04]\n"
     ]
    }
   ],
   "source": [
    "optimizer = ms.nn.Adagrad(molecule_pqcnet.trainable_params(), learning_rate=4e-2)\n",
    "train_pqcnet = ms.nn.TrainOneStepCell(molecule_pqcnet, optimizer)\n",
    "\n",
    "print(\"eps: \", eps)\n",
    "energy_diff = eps * 1000\n",
    "energy_last = initial_energy.asnumpy() + energy_diff\n",
    "iter_idx = 0\n",
    "while abs(energy_diff) > eps:\n",
    "    energy_i = train_pqcnet().asnumpy()\n",
    "    if iter_idx % 5 == 0:\n",
    "        print(\"Step %3d energy %20.16f\" % (iter_idx, float(energy_i)))\n",
    "    energy_diff = energy_last - energy_i\n",
    "    energy_last = energy_i\n",
    "    iter_idx += 1\n",
    "\n",
    "print(\"Optimization completed at step %3d\" % (iter_idx - 1))\n",
    "print(\"Optimized energy: %20.16f\" % (energy_i))\n",
    "print(\"Optimized amplitudes: \\n\", molecule_pqcnet.weight.asnumpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 总结\n",
    "\n",
    "在本案例中，我们通过两种方法，利用量子神经网络得到了LiH分子的基态能量。在第一种方法中，我们利用MindQuantum打包好的`generate_uccsd`函数生成了能够解决该问题的量子神经网络，而在第二种方法中，我们一步一步的构造出了类似的量子神经网络。最终得到的结果是一致的。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "MindSpore",
   "language": "python",
   "name": "mindspore"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}